function [theta11Step3,outputStep3,setupStep3,theta2Step3,AVarTheta2Step3,logDiffSigma] = ...
    step3_theta11_iterate_macro(theta11Step3,theta12OptimalOn,...
    upperBounds,lowerBounds,Insigma,optim,...
    AVarTheta11Step3,outputStep3,bandwidthTstep2,setupStep3,h0RegimeOn,hxRegimeOn,numBootStep2,...
    namesTheta12,theta1Step1,AVarTheta1Step1,setupStep1,theta2Step2)


[theta2Step3,AVarTheta2Step3] = step2SR_threshold_macro(AVarTheta11Step3,outputStep3,bandwidthTstep2,setupStep3.econAct,h0RegimeOn,hxRegimeOn,numBootStep2);
[~,~,~,~,sigmaStep2] = unfoldTheta2_threshold_macro(theta2Step2,setupStep1.nm*2+setupStep1.nn);
[~,~,~,~,sigmaStep3] = unfoldTheta2_threshold_macro(theta2Step3,setupStep3.nm*2+setupStep3.nn);
diff_sigmaStep3      = max(abs(sigmaStep2(:)-sigmaStep3(:)))
logDiffSigma         = nan(5,1);
logDiffSigma(1,1)    = diff_sigmaStep3;   
epsValue = 1D-9;
optim.numOptim = 1;
idx = 0;
while diff_sigmaStep3 > 0.0001 && idx < 4
    idx = idx + 1;
    
    % Re-estimating theta11
    [theta12Step3,AVarTheta12Step3] = step3SR_theta12(namesTheta12,theta1Step1,AVarTheta1Step1,...
        theta2Step3,AVarTheta2Step3,theta12OptimalOn);
    theta11Step3start = theta11Step3;
    for i=1:length(namesTheta12)
        theta11Step3start.(namesTheta12{i}) = 0.01; %not used
    end    
    [~,theta11Step3,outputStep3,setupStep3] = step3SR_theta11(theta12Step3.Robust,...
        namesTheta12,theta11Step3start,...
        setupStep1,upperBounds,lowerBounds,Insigma,optim);

    % Standard errors of theta11: The asymptotic standard errors
    AVarTheta11Step3 = getAVarTheta11Step3(namesTheta12,theta11Step3,AVarTheta1Step1,AVarTheta2Step3,AVarTheta12Step3,...
        setupStep1,setupStep3,epsValue,theta12OptimalOn);

    % updating sigma
    [theta2Step3New,AVarTheta2Step3] = step2SR_threshold_macro(AVarTheta11Step3,...
        outputStep3,bandwidthTstep2,setupStep3.econAct,h0RegimeOn,hxRegimeOn,numBootStep2);
    [~,~,~,~,sigmaStep3]    = unfoldTheta2_threshold_macro(theta2Step3,setupStep3.nm*2+setupStep3.nn);
    [~,~,~,~,sigmaStep3New] = unfoldTheta2_threshold_macro(theta2Step3New,setupStep3.nm*2+setupStep3.nn);
    diff_sigmaStep3 = max(abs(sigmaStep3(:)-sigmaStep3New(:)))
    logDiffSigma(idx+1,1) = diff_sigmaStep3
    theta2Step3 = theta2Step3New;
    
end
